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Why does preferential diffusion strongly 
affect premixed turbulent combustion? 

By V. R. Kuznetsov 1 


1. Introduction 

^ Combustion of premixed reactants in a turbulent flow is a classical but unresolved 
problem. The key problem is to explain the following data: the maximal turbulent 
and laminar burning velocities Uj and ui occur at different equivalence ratios $ . 
(for a review of experimental data, see Kuznetsov & Sabel’nikov (1990), Chapter 
6). Some examples of fuel behavior are: H 2 displays a large shift of the maximum 
value of u t towards the lean mixture, CH 4 has a small lean shift, C 2 H 6 has no shift, 
C 2 H$ has a moderate rich shift, and benzene has a large rich shift. The shift can 
be quite large. For example, the maximum of uj occurs at $ = 10 for H 2 and at 
$ = 1.4 for benzene, while, the maximum of ui is at $ = 1.7 for H 2 and $ = 1 
for benzene. This shift is observed over a large range of Damkohler number, but is 
more pronounced at low Damkohler number. A theory should explain these data. 

It can be seen that the fuels in the above-mentioned sequence are arranged ac- 
cording to the ratio of the molecular diffusivities of oxygen (£>„) and fuel (£>/)• It 
can, therefore, be hypothesized that preferential diffusion strongly affects turbulent 
combustion in all regimes. The correlation between u t and blow-off velocity, based 
on this assumption, is very good over a wide range of conditions (Kuznetsov & 
Sabel’nikov (1990), Chapter 6). 

If the reaction zone were distributed, the influence of molecular diffusivity vari- 
ation should be unimportant since only large-scale fluctuations should affect the 
reaction and their properties do not depend on Reynolds number. On the other 
hand, if the flame front were thin (which was verified by direct numerical simula- 
tion (Rutland et al„ 1989; Trouve, 1991)), the Reynolds number (based on the flame 
front thickness) would be small and the influence of preferential diffusion could be 
significant. 

It is known that the equivalence ratio varies along a curved flame if Df + D 0 . 
However, the mean flame radius of curvature is much larger than the laminar flame 
thickness S L . Therefore, significant influence of preferential diffusion should occur 
only if the flame propagation speed varies with flame curvature. This conclusion 
agrees with Zel’dovich’s long-standing idea about the important role of leading 
points of a flame (the points L U L 2 ,.. . in Fig. 1 which are deep inside the fresh 

mixture). 

The main objective of this paper is to prove Zel’dovich’s hypothesis. An equation 
for the mean flame surface area density (MFSAD) will be employed for this purpose^ 
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FIGURE 1. A sketch of turbulent flame front. Li,L 2 , . . . are leading points. 

a popular treatment (for a review see Candel, et al, 1990). An equation for MFSAD 
can be written 


3E 


dt 


+ <u k > 
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( 1 ) 


where D t is a turbulent diffusivity, H and G are positive functions depending on 
turbulence properties, and E/ is MFSAD. 

- The second objective of this paper is to suggest a different approach to the deriva- 
tion of the equation for MFSAD. It is based on the pdf equation for the reaction 
progress variable C and the relation between the pdf and MFSAD (Kuznetsov & 
Sabel’nikov (1990)). As will be seen later, this treatment suggests an entirely dif- 
ferent closure assumption. 

2. Assumptions 

We wish to prove the hypothesis about the crucial role of leading points for the 
case of equal diffusivities, so we consider the case in which the diffusivities of all 
species and heat are equal. 

The main properties of turbulent flames can be correlated using non-dimensional 
numbers based on two characteristics of the laminar flame, u L and S L (Kuznetsov 
& Sabel’nikov (1990), Chapter 6). If this is the case, the detailed chemical source 
term is irrelevant, and it is important only to model ui and 8i correctly. This 
reasoning is widely used (see Rutland, et al, 1989). It can then be assumed that 
the source term in the equation for the reaction progress variable, C, depends only 
on C itself so that the governing equation is 


dC dc 
p ^ + pUk dT k 


= VDpVC -I- W{C) 


( 2 ) 


We shall assume that the dependence of D or p on C is known. 
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If only these two quantities (i.e., ul and 8l) are important, it is sufficient to 
study the high activation energy limit, i.e., it is assumed that 

W = 0 if C <C a . (3) 


I wdC = const 

Jc 0 


(4) 


where C 0 —* 1. It is clear from Eq. (4) that ui remains constant as Co 1. ^ 

If c o _> the chemical reaction zone is located on the surface on which C = 1 . 
This definition is required since a large volume is occupied by fully burnt products 
in which C = 1. If C 0 is close to unity, the usual analysis of the solution within the 
reaction zone proves that 

oC _ PqUl ( 5 ) 

dn pi,Dk 


Here n is normal to the surface C = 1", and subscripts 0 and b denote fresh mixture 
and burnt products. It should be kept in mind that Eq. (5) is exact even when the 

flamelet model is not valid. ^ . . 

There may be some confusion if the limit C 0 — * 1 is studied. Indeed, Eq. (4) 
shows that W tends to const 6{C- 1). At the same time, there is a term proportional 
to 8(C - 1) in the pdf (because there is a non- zero volume where C = 1). Since 
there is a term d{WP)dC in the equation for the pdf, there can be some confusion. 
To avoid it, let us choose W{C)with{C 0 + 1) so that value C = 1 occurs far from 
the flame. For example, we might let W = const if C 0 < C < 1. For the plane, 
steady laminar flame, one can prove that C = 1 at a finite distance from the flame. 
The equation for the pdf can be obtained, and only then will the limit C 0 — ► 1 be 
investigated. 

It will also be assumed that 8 L is much smaller than the Kolmogorov scale. From 
the mathematical point of view, this means that 6 L -> 0 (or equivalently D -> 0 
since 8l ~ D/ul ) while other quantities are kept constant. 


3. Some geometrical considerations 

According to the Zel’dovich model, reaction takes place only on a surface. Hence 
there are combustion products ( C = 1) on just one side of this surface, and the vol- 
ume of products is not zero, i.e. the pdf of C contains a singular term proportional 
to 8(C - 1), in agreement with the Bray-Moss (1977) model. 

There is another singular term in the Bray-Moss model which is proportional to 
6(C) and represents pockets of fresh mixture. However, this term should be absent 
at any finite Reynolds number (Kuznetsov & Sabel’nikov (1990), Chapter 1) since 
heat spreads simultaneously throughout the whole volume of fluid due to molecular 
diffusion. This means that pockets of fresh mixture are diluted by the products. 
The degree of such a dilution is small at large Reynolds number, i.e. this term arises 
only at infinite Reynolds number. In this limiting case, the flame front should be 
treated as infinitely thin and the continuous term in the Bray-Moss model should 
be absent. It will be assumed in the present treatment that the Reynolds number 
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is finite and, hence, the second singular term is absent. This means that there is 
a minimal value C = C mi „ of the reaction progress variable (which depends on 
distance C m i n from the flame). 

Hence, one has 

p = yS(C — l) + F (6) 


where F is continuous if C min < C < 1 and 7 is the probability of observing the 
value C — 1. Clearly, 7 can be treated as a combustion efficiency if the flamelet 
model is valid. 

To calculate the continuous term in Eq. (6), one can use the equation (Cant et 
al , 1990) 


d(Yi\C — cons t) _ .dC 
dv ' dn 


cons 


*) 


(7) 


where E is the area of surface on which C — const (so the left-hand side of Eq. (7) 
is the mean surface area density). Surfaces with different valuts of C have almost 
the same area since the thickness of the flame is small compared to the Kolmogorov 
scale. Hence, 


d(ll\C = const ) 
dv 


£/(f). 


( 8 ) 


To calculate the conditionally averaged gradient of the reaction progress variable, 
let us choose some point on the surface C — 1” and the frame moving with the 
velocity of this point. Then, at a distance much smaller than the Kolmogorov scale 
and much larger than <5^, one sees an almost plane, steady laminar flame, so that 
Eq. (2) reduces to 


PoUL 


dC_ 

dn 


d dC 

= Tn D ^ + W 


which is familiar from laminar flame theory. If C„ — ► 1, the chemical source term 
can be neglected. Hence, after the integration one has 


Po^lC = Dp 


dC 

dn 


There are no random parameters, thus 


dC 

dn 


C = const ) = — C 


Dp 


(9) 


( 10 ) 


Combining Eqs. (7), (8), and (10), one has 


F _ y 

p 0 u L C f 


( 11 ) 


It should be kept in mind that Eq. (11) is exact for C = 1~ since Eq. (5) is exact 
if C 0 1. Eq. (11) is approximate for C < 1 since any influence of turbulence on 
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flame structure was neglected in Eq. (9). This influence is not necessarily small for 
small C since surfaces C = const ~ C min are located far from the flame where the 
structure of the scalar field is heavily affected by the turbulence. This conclusion 
agrees with direct numerical simulations (Rutland et al, 1989) performed at large 
Damkohler number. 

This is not a serious deficiency of Eq. (11) since the mean value of any function 
if(C) can be easily calculated using identities 

E/ P rk Y>( 

< <p >=< 0) > +¥>(0) = 7V(1)+(1-tM0)+— J o C dC 

The integrand has no singular points, i.e. the domain C ~ 0 where Eq. (11) is 
not valid does not play a significant role. 

One can see that the pdf of reaction progress variable depends only on two func- 
tions of coordinates, the mean flame surface area density E / and the combustion 
efficiency 7 . 


4. Exact equations for E/ and 7 

It is natural to try to obtain equations for E / and 7 using an exact (but unclosed) 
equation for the pdf. Using methods developed by Kuznetsov & Sabel’nikov (1990), 
one can obtain two equivalent forms of the pdf equation. 

£-1 pv k (Z , C)F + pv k (x, 1)7 S(C - 1)] = --^N(x,C)pF - -~WF (12) 

-^-[pv k (x,C)F + pv k (S, 1)^(0 - 1)] = --^FA(x,C) - ^WF (13) 

where v k is a flow velocity averaged at the condition C = const, N and A are 
quantities D(dC/dx k ) 2 and VDpVC averaged at the same condition. These equa- 
tions are valid only for the chemistry model adopted in Section 2. The presence 
of a ^-function on the left-hand side can lead to some confusion. To clarify it, let 
us note the F = 0 if C > 1 and F # 0 if C < 1. It can be guessed that another 
^-function will appear on the right-hand side on differentiation of the function F. 

Quantities N and A can be calculated to first approximation using the flamelet 
model. For example, using Eq. (9), one has 


N = 


pWl . 

p 2 D 


(14) 


Note that this equation is exact 
Similarly, one has 


if C = 1 and C Q 


A== 

Dp 


1 since Eq. (5) is exact. 


(15) 
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The flamelet model is valid if D -*■ 0 and other quantities are kept constant. It 
is seen from Eq. (15) that A — ► oc if D — > 0. Therefore, Eq. (13) reduces to 

^-FA = 0 if C<C 0 . 

OC 

This equation is easily integrated. After the use of Eq. (15), one again has Eq. 
(11), which was obtained previously using purely geometrical considerations. 

Thus Eq. (11) is the first term of the asymptotic expansion of the solution of 
the pdf equation. The difficulty is that to calculate MFSAD (which appears as a 
constant of integration), other terms must be investigated, i.e. small, previously 
neglected terms in Eq. (14) or (15) should be taken into account. This means that 
the small influence of turbulence on the inner structure of the flame front, flame 
stretch, and flame curvature should be included. Clearly, in the general case, it will 
also include effects of preferential diffusion. 

There is a simpler approach based on integral relations which can be deduced 
from the exact Eqs. (12) and (13) without use of the approximate relations (14) 
and (15). In performing the calculations, one should keep in mind that the pdf is 
defined for all values of C( — oo < C < oo) and 


F = 0 if C<0 or C>1. (16) 

To obtain the first equation, let us integrate Eq. (13) from C 0 to C = oo. Using 
Eqs. (3) and (16) and assuming that C a — ► 0, one has 

Q 

Q^-pbV,{x, 1)7 = F(x, l")A(f, 1). 

Small corrections to Eq. (15) do not play a significant role here. Hence, using Eqs. 
(11) and (15), one has 

d ^ , 

— p b Vi(x, 1)7 = PoUlT, f. (17) 

To obtain the second equation, one should multiply Eq. (12) by C and integrate 
from C = -oo to C = C 0 . Using Eqs. (3) and (16), one has 

jf ' pv,FdC = [NFp - C-^NFp] c =c 0 . 


Assuming C 0 — ♦ 1 and using Eq. (14), which is exact for C = 1, one has 


d_ 

dr, 



pv,CFdC 


pWl 

[ Dbpb 


F- 


dNpF' 

dC 


c= 1 


(18) 


Note that the second term on the right-hand side of Eq. (18) can not be calculated 
exactly using Eqs. (11) and (14). To calculate it, let us compare Eqs. (12) and 
(13). One has 
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Integrating, one has 

^N„F = FA. (19) 

The term FA can not be calculated exactly from Eqs. (11) and (15). However, 
the use of Eq. (19) has some advantages since A is related to a mean velocity of 
chemical reaction zone which has a clear physical meaning. Since the reaction zone 
is located where C = 1", its velocity in a laboratory frame u\ is given by 


u 


ch 

it — 


dc dc / (dc\ 

dt dxk / 


if C = Co- 


Using Eqs. (2) and (20), one has 


u 


ch __ 
k — 


A dC /(dC_ \ 2 

Pkdx k / 


- Uk 


if 


C = C 0 . 


( 20 ) 


( 21 ) 


Evidently, Eq. (21) is valid if C 0 -> 1. So the normal component of chemical 
reaction zone velocity relative to the local velocity of burnt products is 



It should be kept in mind that u> is negative (if the normal to the flame is directed 
toward the fresh mixture). To avoid confusion, it is convenient to change the sign 
in Eq. (22) so that u is positive. Using Eq. (15), one has 


= -^-A. (23) 

poUL 


This equation is valid on the surface C = 1" for instantaneous values of u and A. 
At the same time, one needs to calculate the conditionally averaged value of A. This 
means that in the homogeneous case, quantity A should be averaged over a volume 
between surfaces C = 1“ and C = 1 - dC. Generally speaking, the averages over 
this volume and over the surface C = 1" are different due to the random fluctuation 
of the distance between the surfaces C = 1~ and C = 1 - dC. However, in the 
studied case, such fluctuations are absent since this distance is 


“ (dC/du) 

where dC/du is constant according to Eq. (5). Hence, Eq. (22) can be averaged 
over surface C = 1" to calculate the conditionally averaged value of A. 

Using Eqs. (18), (19), and (23), one has 


d_ 

dx, 


i: 


pviCFdC 


7^(/>o«l - /»tw)F|c=i- 

Vbpb 
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The approximate Eq. (11) can now be used to obtain 


d 

dx, 



VidCT,f = Po«l(0U£, — w)E/ 


( 24 ) 


where 0 = p 0 /pb- 

This equation is exact if the flamelet model is used. It is seen that two quantities 
should be modeled: the conditionally averaged gas velocity and the difference Oul — 
lo. It is clear that Ou/, is the velocity of a plane steady laminar flame relative to 
the burnt products. Hence, u would equal to 0u L if the properties of a flame front 
in a turbulent flow were exactly the same as those of a laminar flame. Of course, 
this is not true. In fact, a small difference between 0ujr, and w is caused by the 
flame curvature and strain, which can not be neglected since the left-hand side of 
Eq. (23) is also small (note that D — ► 0). 

5. Turbulent diffusion in a flame 

Consider now the quantity v k . It is convenient to begin with the constant density 
case. Since the flame front is very thin, there is no appreciable velocity difference 
across it, i.e. Vk does not depend on C. In a variable density case, the variation of 
gas velocity across the flame front can be calculated from the continuity equation. 
Since the flame front is almost planar, in a frame moving with the front 

P 

Uk = — u L n k 

Po 

where n k is a unit vector normal to a flame front. Hence, in the laboratory frame, 

w* = Uk(£, l - ) + —UlUk 

Po 

where u k {x , l - ) is the gas velocity on the surface C = 1 ~. After averaging at the 
condition C = 1", one has 


v k = Vk{x, 1 ) + - — —ui<rik\C=\ >. 
Po 


(25) 


Consider now the case u' >> u /, where u' is the rms flow velocity. Then the 
second term on the right-hand side of Eq. (25) can be neglected. Also, the flame 
moves with almost the same velocity as the gas. Hence, the theory of turbulent 
diffusion can be applied. This means that the quantity v k should be treated as 
differential operator 


Vk(x, 1 ) =< Uk > —D, 


d 

dx k 


( 26 ) 


where D t is a turbulent diffusivity. Using Eqs. (24) and (26), one has 


9 < Uk > £/ _ d n d£/ Po(Qm. ~ u))u L ^ 
dx k dx k dx k Dp 2 dC ~‘ f 


( 27 ) 


451 


Preferential diffusion effects on premixed turbulent combustion 

It is worth noting that the validity of Eq. (27) is questionable when u' ~ u L 
because the surface C = 1“ is moving rapidly relative to the gas. At the same time, 
the surface on which a non-reactive scalar is constant moves with the Kolmogorov 
velocity relative to the gas. This is easily seen from Eq. (22) which is also valid for 
a non-reacting scalar if u c k h is the velocity of an isoscalar surface. Obviously, the 
quantities A and ( dC/dxif have Kolmogorov scaling. Hence, the left-hand side of 
Eq. (22) also has Kolmogorov scaling, i.e. ui « u' . 

Therefore, it could be expected that if u' ~ u L , the quantity D t in Eq. (27) is 
not the turbulent diffusivity of non-reactive scalar. Even counter-gradient diffusion 
can be expected. 

The same approach can be used to close Eq. (17)> be. 


d 9 n ^7 a v 


(28) 


6. Effects of flame curvature 

As was mentioned earlier, one can not neglect the difference between Out and 
w. It is easily seen that if one assumes Out - w = const - const E/, then Eq. (27) 
reduces to Eq. (1). However, we do not intend to develop a model for this difference 
since, as will be seen later, we can reach important conclusions by analyzing only 
the sign of difference Out - w in various parts of a flame. 

We shall consider only the effects of flame curvature. Obviously, flame curvature 
in low and high mean temperature regions have opposite predominant signs. This 
becomes clear from examination of Fig. 1 where a flame homogeneous in 12 direction 
is shown. In a low mean temperature region (x 2 < 0), the flame front xi(x 2 ) has 
more minima than maxima. In a high mean temperature region (x 2 > 0), the flame 
front has more maxima than minima. One more example is the Gaussian curve 
^i(x 2 ). In this case, one has 


( 


<Px\ 
dx 2 


Xi 


= const 


) 



*i 

<x> >’ 


i.e. the flame curvature is positive in a low mean temperature region and is negative 
in a high mean temperature region. 

It was shown long ago (Markstein, 1964) that ©u L > w if the curvature is negative 
and Qui < uj if the curvature is positive. So Out < in a high mean temperature 
part of a flame and 0 u L > u> in a low mean temperature part of a flame. For the 
plane turbulent flame shown in Fig. 1, one has 

0u/, — u> > 0 if xi — + —00 (28) 


0U£ — u> < 0 if xj — » + 00 . 


(29) 
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7. Plane turbulent flame 

Consider now a steady plane turbulent flame such as the one shown in Fig. 1. 
One has 

< u\ > = ut if x\ —* — oo 

< u\ > = Qut if xj — ► +oo. 

Let us analyze solutions of Eq. (27) in the regions xi — » +oo and xi — ► -oo, 
assuming that all coefficients in this equation tend to constants in these limits. 
One has 


£/ = const ex p(^" (©Ut ± ~ ^ | 


if 


X! 


+oo 


where 

_ 4DiP^(Qul — u)ul 

lo Dp 2 dc 

One can see from Eq. (29) that w < 0; hence, one solution increases without 
bound and the other tends to zero as Z] — * oo. 

One has also 

E/ = const exp| (u ( ± \Ju 2 - w 

One can see from Eq. (28) that w > 0. So there are two cases. The first one is 
Ut < y/w, in which case each solution changes sign and is meaningless. The second 
case is u t > y/w, and both solutions tend to zero as ii -» — oo. Therefore, the 
spectrum of eigenvalues is continuous. Such spectra are often met in the combustion 
theory and other areas (Barenblatt, et al, 1985). It was proven that the physical 
meaningful velocity is the minimal value of flame velocity (the other values are due 
to the initial disturbances). Hence, one has 



u t 


( 4Z>u/,p^(0U£, — w) 

To Dp 2 dC 


Since u t depends only on value w deep in the fresh mixture, this proves Zel’dovich’s 
assumption about the crucial role of leading points. 

As was said above, the difference 0U£ - u> depends on flame curvature. In the 
general case, the value of this difference depends on effects caused by the preferential 
diffusion. This means that 0 u — w and u l should be calculated at an equivalence 
ratio different from the overall equivalence ratio. This is the cause of a shift between 
maximal values of u t and ui- 

The sign of this shift can be established by studying simple configuration (say, 
a spherical flame). It was found that < $ if D t > D a and $/ > $ if Df < 
D 0 (Kuznetsov & Sabel’nikov, 1990). Here, $ f is the equivalence ratio in the 
reaction zone and $ is the overall equivalence ratio. Therefore, a lean shift should 
be observed in # 2 -air flame and a rich shift should be observed in C 3 H & - air flame. 
This conclusion is fully consistent with experimental data (see Introduction). 
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8. Conclusion 

This investigation led to three important conclusions. 

i) The influence of turbulence on the inner structure of a flame front is quite impor- 
tant even when the laminar flame thickness is much smaller than the Kolmogorov 

ii) Turbulent flame propagation depends on conditions within the low temperature 
part of a flame. Preferential diffusion plays an important role here because flame 
curvature in this part of the flame has a predominant sign. 

iii) To develop a correct equation for the mean flame surface area density, one should 
properly model the difference between laminar flame velocity and the actual ve- 
locity of a flame front. 
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